Strong-coupling expansion for the two-species Bose-Hubbard model 
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To analyze the ground-state phase diagram of Bose-Bose mixtures loaded into d-dimensional hypercubic op- 
tical lattices, we perform a strong-coupling power-series expansion in the kinetic energy term (plus a scaling 
analysis) for the two-species Bose-Hubbard model with onsite boson-boson interactions. We consider both 
repulsive and attractive interspecies interaction, and obtain an analytical expression for the phase boundary be- 
tween the incompressible Mott insulator and the compressible superfluid phase up to third order in the hoppings. 
In particular, we find a re-entrant quantum phase transition from paired superfluid (superfluidity of composite 
bosons, i.e. Bose-Bose pairs) to Mott insulator and again to a paired superfluid in all one, two and three di- 
mensions, when the interspecies interaction is sufficiently large and attractive. We hope that some of our results 
could be tested with ultracold atomic systems. 

PACS numbers: 03.75.-b, 37.10.Jk, 67.85.-d 



I. INTRODUCTION 

Single-species Bose-Hubbard (BH) model is the bosonic 
generalization of the Hubbard model, and was introduced 
originally to describe "^He in porous media or disordered gran- 
ular superconductors [1]. For hypercubic lattices in all di- 
mensions d, there are only two phases in this model: an in- 
compressible Mott insulator at commensurate (integer) fill- 
ings and a compressible superfluid phase otherwise. The su- 
perfluid phase is well described by weak-coupling theories, 
but the insulating phase is a strong-coupling phenomenon that 
only appears when the system is on a lattice. Transition from 
the Mott insulator to the superfluid phase occurs as the hop- 
ping, particle-particle interaction, or the chemical potential is 
varied fT]. 

It is the recent observation of this transition in effectively 
three- [2], one- [3], and two-dimensional 10, HI] optical lat- 
tices, which has been considered one of the most remarkable 
achievements in the field of ultracold atomic gases, since it 
paved the way for studying other strongly correlated phases 
in similar setups. Such lattices are created by the intersection 
of laser fields, and they are nondissipative periodic potential 
energy surfaces for the atoms. Motivated by this success in 
experimentally simulating the single-species BH model with 
ultracold atomic Bose gases loaded into optical lattices, there 
has been recently an intense theoretical activity in analyzing 
BH as well as Fermi-Hubbard type models |6]. 

For instance, in addition to the Mott insulator and single- 
species superfluid phases, it has been predicted that the two- 
species BH model has at least two additional phases: an in- 
compressible super-counter flow and a compressible paired 
superfluid phase ff^T^. Our main interest here is in the latter 
phase, where a direct transition from the Mott insulator to the 
paired superfluid phase (superfluidity of composite bosons, 
i.e. Bose-Bose pairs) has been predicted, when both species 
have integer fillings and the interspecies interaction is suffi- 
ciently large and attractive. Given that the interspecies inter- 
actions can be fine tuned in ongoing experiments, e.g. with 
4iK-87Rb with mixtures iflTi [Tsll . via using Feshbach reso- 
nances, we hope that some of our results could be tested with 
ultracold atomic systems. 



In this paper, we examine the ground-state phase diagram of 
the two-species BH model with on-site boson-boson interac- 
tions in d-dimensional hypercubic lattices, including both the 
repulsive and attractive interspecies interaction, via a strong- 
coupling perturbation theory in the hopping. We carry the 
expansion out to third-order in the hopping, and perform a 
scaling analysis using the known critical behavior at the tip 
of the insulating lobes, which allows us to accurately predict 
the critical point, and the shape of the insulating lobes in the 
plane of the chemical potential and the hopping. This tech- 
nique was previously used to discuss the phase diagram of the 
single-species BH model lfT9l - l23ll . extended BH model ll24ll . 
and of the hardcore BH model with a superlattice llzsll . and 
its results showed an excellent agreement with Monte Carlo 
simulations ll23l Izsll . Motivated by the success of this tech- 
nique with these models, here we apply it to the two-species 
BH model, hoping to develop an analytical approach which 
could be as accurate as the numerical ones. 

The remaining paper is organized as follows. After in- 
troducing the model Hamiltonian in Sec. Ull we develop the 
strong-coupling expansion in Sec. |III1 where we derive an 
analytical expression for the phase boundary between the in- 
compressible Mott insulator and the compressible superfluid 
phase. Then, in Sec. lIVI we propose a chemical-potential ex- 
trapolation technique based on scaling theory to extrapolate 
our third-order power-series expansion into a functional form 
that is appropriate for the Mott lobes, and use it to obtain typ- 
ical ground-state phase diagrams. A brief summary of our 
conclusions is given in Sec.lV] 



II. TWO-SPECIES BOSE-HUBBARD MODEL 

To describe Bose-Bose mixtures loaded into optical lattices, 
we consider the following two-species BH Hamiltonian, 
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where the pseudo-spin a = {f, i} labels the trapped hyper- 
fine states of a given species of bosons, or labels different 
types of bosons in a two-species mixture, t^ cr is the tun- 



A. Ground-State Wave Functions 



The perturbation theory is performed with respect to the 



neling (or hopping) matrix between sites i and j, b\_^ {bi,a) ground state of the system when ^ = 0, and therefore 
is the boson creation (annihilation) and ni,„ = 6| _6i.cr i: 



the boson number operator at site i, Uaa' is the strength of 
the onsite boson-boson interaction between a and a' compo- 
nents, and /icr is the chemical potential. In this manuscript, 
we consider a d-dimensional hypercubic lattice with M sites, 
for which we assume t^j ^ is a real symmetric matrix with el- 
ements cr = icr > for i and j nearest neighbors and 
otherwise. The lattice coordination number (or the number of 
nearest neighbors) for such lattices is z = 2d. 

We take the intraspecies interactions to be repulsive 
({C/ff , Ui^} > 0), but discuss both repulsive and attractive 
interspecies interaction U^i as long as U^^Un > U^^. This 
guarantees the stability of the mixture against collapse when 
U^i ^ 0, and against phase separation when U^i 0. How- 
ever, when the interspecies interaction is sufficiently large and 
attractive, we note that instead of a direct transition from the 
Mott insulator to a single particle superfluid phase, it is possi- 
ble to have a transition from the Mott insulator to a paired su- 
perfluid phase (superfluidity of composite bosons, i.e. Bose- 
Bose pairs) jT Hlql . Therefore, one needs to consider both 
possibilities, as discussed next. 



I 

we first need zeroth order wave functions of the Mott phase 
and of its defect states. To zeroth order in and t^, the Mott 
insulator wave function can be written as. 
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where (rii^a-) = is an integer number corresponding to the 
ground-state occupancy of the pseudo-spin a bosons, (• • • ) is 
the thermal average, and |0) is the vacuum state. On the other 
hand, the wave functions of the defect states are determined 
by degenerate perturbation theory. The reason for that lies 
in the fact that when exactly one extra elementary particle or 
hole is added to the Mott phase, it could go to any of the M 
lattice sites, since all of those states share the same energy 
when = tj^ = 0. Therefore, the initial degeneracy of the 
defect states is of order M. 

When the elementary excitations involve a single-a-particle 
(exactly one extra pseudo-spin a boson) or a single-a-hole 
(exactly one less pseudo-spin a boson), this degeneracy is 
lifted at first order in and ij^. The treatment for this case is 
very similar to the single-species BH model lfT9il24.1 , and the 
wave functions (to zeroth order in and tj^) for the single-a- 
particle and single-cr-hole defect states turn out to be 



III. STRONG-COUPLING EXPANSION 
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We use the many-body version of Rayleigh-Schrodinger 
perturbation theory in the kinetic energy term to perform the 
expansion (in powers of and t]) for the different energies 
needed to carry out our analysis. The strong-coupling expan- 
sion technique was previously used to discuss the phase di- 
agram of the single-species BH model lll9l - l2lll23ll . extended 
BH model i24tl . and of the hardcore BH model with a superlat- 
tice [iH, and its results showed an excellent agreement with 
Monte Carlo simulations ll23l Esll . Motivated by the success 
of this technique with these models, here we apply it to the 
two-species BH model. 

To determine the phase boundary separating the incom- 
pressible Mott phase from the compressible superfluid phase 
within the strong-coupling expansion method, one needs the 
energy of the Mott phase and of its 'defect' states (those states 
which have exactly one extra elementary particle or hole about 
the ground state) as a function of and t^. At the point 
where the energy of the incompressible state becomes equal 
to its defect state, the system becomes compressible, assum- 
ing that the compressibility approaches zero continuously at 
the phase boundary. Here, we remark that this technique can- 
not be used to calculate the phase boundary between two com- 
pressible phases. 



where 



jsah j.jjg eigenvector of the hopping matrix 



tij.cr with the highest eigenvalue (which is zta with z = 2d) 
such that J2j tij.a-fj'^^ = ztcrfl'^^ ■ The normalization condi- 
tion requires that l/f^^P — 1. Notice that we choose the 
highest eigenvalue of tij,„ because the hopping matrix enters 
the Hamiltonian as —tij_a, and we ultimately want the lowest- 
energy states. 

However, when the elementary excitations involve two par- 
ticles (exactly one extra boson of each species) or two holes 
(exactly one less boson of each species), the degeneracy is 
lifted at second order in and i^. Such elementary excita- 
tions occur when U-\i is sufficiently large and attractive [2^, 
and the wave functions (to zeroth order in and t]) for the 
two-particle and two-hole defect states can be written as 
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/* turns out to be the eigenvector of the 



tij,t^ij,i matrix with the highest eigenvalue (which is zt^t^ 
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Since the 



with z = 2d) such that J2j iy\t^y,i// = ^^t^i/, 
elementary excitations involve two particles or two holes, the 
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degenerate defect states cannot be connected by one hopping, 
but rather requke two hoppings to be connected. Therefore, 
one expects the degeneracy to be Hfted at least at second order 
in and ij., as discussed next. 



B. Ground-State Energies 

Next, we employ the many-body version of Rayleigh- 
Schrodinger perturbation theory in t-^ and with respect to 
the ground state of the system when = ti = 0, and cal- 
culate the energy of the Mott phase and of its defect states. 
The energy of the Mott state is obtained via nondegenerate 
perturbation theory, and to third order in and it is given 
by 



M 



(7) 



This is an extensive quantity, i.e. iS^ott proportional to the 
number of lattice sites A/. The odd-order terms in and 
vanish for the d-dimensional hypercubic lattices considered in 
this manuscript, which is simply because the Mott state given 
in Eq. (|2]i cannot be connected to itself by only one hopping, 
but rather requires two hoppings to be connected. Notice that 
Eq. O recovers the known result for the single-species BH 
model when one of the pseudo-spin components have vanish- 
ing filling, e.g. 71^ = ESIH. 



The calculation of the defect-state energies is more involved 
since it requires using degenerate perturbation theory. As 
mentioned above, when the elementary excitations involve a 
single-CT-particle or a single-cr-hole, the degeneracy is lifted at 
first order in and t^. A lengthy but straightforward calcula- 
tion leads to the energy of the single-cr-particle defect state up 
to third order in and as 
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where (— f) =1 and vice versa. Here, we assume C/o-o- ^ ^o- and {TJ-a-a, ± U^l\} ^ ^-o-- Equation (|8]l is valid for 

all d-dimensional hypercubic lattices, and it recovers the known result for the single species BH model when ri-a- = Li9„24il. 
Note that this expression also recovers the known result for the single species BH model when TJ^^ = 0, which provides an 
independent check of the algebra. To third order in and t^, we obtain a similar expression for the energy of the single-cr-hole 
defect state given by 
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which is also valid for all d-dimensional hypercubic lattices, and it also recovers the known result for the single-species BH 
model when n^a- = or U^^ = 1 19.,24il . Here, we again assume TIa-a- ^ and {TJ-a-a, \TJ-a-a ± U^\\] ^ <_cr- We also 
checked the accuracy of the second-order terms in Eqs. ([Hi and (|9]) via exact small-cluster (two-site) calculations with one a and 
two — cr particles. 

We note that the mean-field phase boundary between the Mott phase and its single-cr-particle and single-cr-hole defect states 
can be calculated as 



par, hoi 



= C/o^K - 1/2) + C/t;n_, - zt,/2 ± ^UIJA - + l/2)zt, + ^2^2/4. 
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This expression is exact for infinite-dimensional hypercubic lattices, and it recovers the known result for the single species BH 
model when n^a = or U^i = 1 1]. In the d — > oo limit (while keeping dt„ constant), we checked that our strong-coupling 
perturbation results given in Eqs. (O and (|9]l agree with this exact solution when the latter is expanded out to third order in and 
t^, providing an independent check of the algebra. Equation dTol ) also shows that, for infinite-dimensional lattices, the Mott lobes 
are separated by U^in^cr, but their shapes and critical points (the latter are obtained by setting = Mct°') are independent of 
[/-[■J,. This is not the case for finite-dimensional lattices as can be clearly seen from our results. It is also important to mention 
here that both the shapes and critical points are independent of the sign of U^i in finite dimensions (at the third-order presented 
here) as can be seen in Eqs. ([Sj and (|9]l. 

However, when the elementary excitations involve two particles or two holes (which occurs when U^i is sufficiently large 
and attractive |.26i1 ). the degeneracy is lifted at second order in and t^. A lengthy but straightforward calculation leads to the 
energy of the two-particle defect state up to third order in and t as 



^dcf = ^M^tt + C^ti("t + + 1) + E (t^-'^- - 

(T 
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Here, we assume {Uaa, \U^i\, 2J7o-ct + U^i} ^ ^<t- Equation ( fTTT i is valid for all d-dimensional hypercubic lattices, where the 
odd-order terms in and vanish l'27'l . To third order in ^ and t^, we obtain a similar expression for the energy of the two-hole 
defect state given by 
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E 



nl _ {nl - 1) 2n^{jia 
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which is also valid for all d-dimensional hypercubic lattices, 
where the odd-order terms in and vanish l27ll . Here, 
we again assume {Uccr,\U^i\,2Uccr + U-\i) > t^. Since 
the single-cr-particle and single-cr-hole defect states have cor- 
rections to first order in the hopping, while the two-particle 
and two-hole defect states have corrections to second order 
in the hopping, the slopes of the Mott lobes are finite as 
{Hi'ti) — >^ in the former case, but they vanish in the lat- 
ter case. Hence, the shape of the insulating lobes are expected 
to be very different for two-particle or two-hole excitations. 
In addition, the chemical-potential widths (/i^r) of all Mott 
lobes are Uaa in the former case, but they {{^^ + Ijli)/2] are 
U^i + {Un + t^u)/2 in the latter. 

We note that in the limit when = = t,U^^ = = 
Uo, = U', = ni = uq, ^ ^ ^ ^, and z ^ 2 
(or d = 1), Eq. (fT2] i is in complete agreement with Eq. (3) of 
Ref. flT'l, providing an independent check of the algebra. In 
addition, in the limit when t-^ — — J, U^^ — = U, 
U-^l = W K -U, Uf = ni = m, and = ALL = A*' 
Eqs. (fTTT i and ( fT2] i reduce to those given in Ref. fl2\ (after 
setting Unn = there). However, the terms that are propor- 
tional to t^ti are not included in their definitions of the two- 
particle and two-hole excitation gaps. We also checked the 
accuracy of Eqs. ( fTTT i and (fT2b via exact small-cluster (two- 
site) calculations with one particle of each species. 

We would also like to remark in passing that the energy dif- 
ference between the Mott phase and its defect states determine 
the phase boundary of the particle and hole branches. This is 
because at the point where the energy of the incompressible 



state becomes equal to its defect state, the system becomes 
compressible, assuming that the compressibility approaches 
zero continuously at the phase boundary. While £'Mott ^^'^ 
its defects E^f, E^^^f, Ef^^ and E^^^^ depend on the lattice 
size M, their difference do not. Therefore, the chemical po- 
tentials that determine the particle and hole branches are inde- 
pendent of AI at the phase boundaries. This indicates that the 
numerical Monte Carlo simulations should not have a strong 
dependence on M. 

It is known that the third-order strong-coupling expansion 
is not very accurate near the tip of the Mott lobes, as and ti 
are not very small there lfl9il24ll . For this reason, an extrapola- 
tion technique is highly desirable to determine more accurate 
phase diagrams. Therefore, having discussed the third-order 
strong-coupling expansion for a general two-species Bose- 
Bose mixtures with arbitary hoppings t^-, interactions Uaa', 
densities n^, and chemical potentials next we show how 
to develop a scahng theory. 



IV. EXTRAPOLATION TECHNIQUE 

In this section, we propose a chemical potential extrapo- 
lation technique based on scaling theory to extrapolate our 
third-order power- series expansion into a functional form that 
is appropriate for the entire Mott lobes. It is known that the 
critical point at the tip of the lobes has the scaling behavior of 
a(d+ l)-dimensional XY model, and therefore the lobes have 
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Kosterlitz-Thouless shapes for d — 1 and power-law shapes 
for d > 1. For illustration purposes, here we analyze only 
the latter case, but this technique can be easily adapted to the 
d = 1 case lfl9ll . 



A. Scaling Ansatz 

From now on we consider a two-species mixture with = 
tl = t, Uf^ = Uii = U, U^i = V, = ni = n, and 
Mt ~ l^i = When d> 1, we propose the following ansatz 
which includes the known power-law critical behavior of the 
tip of the lobes 



(13) 



where 74(a;) = a + hx + cx'^ + dx^ + ■ ■ ■ andi?(x) = a + /3x + 
^x^ + 5x^ + • • • are regular functions of a; = 2dt/U, Xc is the 
critical point which determines the location of the lobes, and 
zi^ is the critical exponent for the {d + l)-dimensional XY 
model which determines the shape of the lobes near Xc = 
2dtc/U. In Eq. (HjI ). the plus sign corresponds to the particle 
branch, and the minus sign corresponds to the hole branch. 
The form of the ansatz is taken to be the same for both single- 
and two-partice (or single- and two-hole) excitations, but the 
parameters are very different. 

The parameters a, b, c and d depend on U, V and n, and 
they are determined by matching them with the coefficients 
given by our third-order expansion such that A{x) = (/iP**' + 
/i'^°*)/(2C/). Here, /^p^'' and are our strong-coupling ex- 
pansion results determined from Eqs. dHJ and (|9]l for the 
single-particle and single-hole excitations, or from Eqs. (fTTT l 
and ( fT2] l for the two-particle and two-hole excitations, respec- 
tively. Writing our strong-coupling expansion results for the 
particle and hole branches in the form /^p**' ~ U J2n=o ^n^" 
and = UJ2n=o^nX"' leads to a = {e^ 

b^{e+ + er)/2, c = (e+ + e^)/2, and d - (e+ 



So")/2, 



n/2. 

To determine the U, V and n dependence of the parameters 
a, /3, 7, 6, Xc and zi', we first expand the left hand side of 
B{x){xc - xY" = (/iP^'' - /i'>°i)/(2C/) in powers of x, and 
match the coefficients with the coefficients given by our third- 
order expansion, leading to 
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We fix ZV at its well-known values such that zv « 2/3 for 
d — 2 and zv — 1/2 for d > 2. If the exact value of Xc 



is known via other means, e.g. numerical simulations, a, /3, 
7 and b can be calculated accordingly, for which the extrap- 
olation technique gives very accurate results |23[ l25il . If the 
exact value of Xc is not known, then we set (5 = 0, and solve 
Eqs. (fl4] i. (fTsT i, ( fTSI l and the 5 = equation to determine 
a, /?, 7 and Xn self-consistently, which also leads to accurate 
results ifioi I24I1 . Next we present typical ground-state phase 
diagrams for [d = 2)- and (d = 3)-dimensional hypercubic 
lattices obtained from this extrapolation technique. 



B. Numerical Results 

In Figs. [U and |2l the results of the third-order strong- 
coupling expansion (dotted lines) are compared to those of the 
extrapolation technique (hollow pink-squares and solid black- 
circles) when V — 0.5U and V = — 0.85[/, respectively, in 
two (d = 2 or z = 4) and three (d = 3 or z = 6) dimensions. 
We recall here that 1^=1^^ t, U^^ = = U, C/tj, = V, 
n^ = n^ = n, and = = ^. 

In Fig. [1] we show the chemical potential fj, (in units of U) 
versus x = 2dt/U phase diagram for (a) two-dimensional and 
(b) three-dimensional hypercubic lattices, where we choose 
the interspecies interaction to be repulsive V = 0.5U. Com- 
paring Eqs. (Is) and (|9]l with Eqs. (fTTT l and (fT2] i. we expect 
that the excited state of the system to be the usual superfluid 
for all > for all t. The dotted lines correspond to phase 
boundary for the Mott insulator to superfluid state as deter- 
mined from the third-order strong-coupling expansion, and 
the hollow pink-squares correspond to the extrapolation fits 
for the single-particle and single-hole excitations discussed in 
the text. We recall here that an incompressible super-counter 
flow phase Ill-Ill [HI] also exists outside of the Mott insulator 
lobes, but our current formalism cannot be used to locate its 
phase boundary. 



TABLE I. List of the critical points (location of the tips) Xc — 
2dtc/U for the first two Mott insulator lobes that are found from 
the chemical potential extrapolation technique described in the text. 
Here, = = t, Uff = Un = U, U^i = V , n-f = ni — n, and 
fj,f = ^il = yL. These critical points for the single-particle or single- 
hole excitations are determined from Eqs. ^ and (|9}, and they tend 
to move in as V increases, and are independent of the sign of V . 
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(a) Two dimensions (V=0.5U) 
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FIG. 1. (Color online) Chemical potential fi (in units of U) versus 
X = 2dt/U phase diagram for (a) two- and (b) three-dimensional 
hypercubic lattices with tf — — t, Uff — C/4.4. = U, Ufi = 
V = 0.5U, n-f = — n, and fif — fi^ = ^. The dotted lines 
correspond to phase boundary for the Mott insulator to superfluid 
state as determined from the third-order strong-coupling expansion, 
and the hollow pink-squares to the extrapolation fit for the single- 
particle or single-hole excitations discussed in the text. Recall that 
an incompressible super-counter flow phase also exists outside of the 
Mott insulator lobes. 



At t = 0, the chemical potential width of all Mott lobes 
are U (similar to the single-species BH model), but they are 
separated from each other by F as a function of As < in- 
creases from zero, the range of /i about which the ground state 
is a Mott insulator decreases, and the Mott insulator phase 



disappears at a critical value of t, beyond which the system 
becomes a superfluid. In additionj^iniilar to what was found 
for the single-species BH model fl9U24tl . the strong-coupling 
expansion overestimates the phase boundaries, and it leads to 
unphysical pointed tips for all Mott lobes, which is expected 
since a finite-order expansion cannot describe the physics of 
the critical point correctly. A short list of V/U versus the crit- 
ical points Xc = 2dtc/U is presented for the first two Mott 
insulator lobes in Table |I] where it is shown that the criti- 
cal points tend to move in as V increases. This is because 
presence of a second species (say —a ones) screens the on- 
site intraspecies repulsion Uaa between cr-species, and hence 
increases the superfluid region. 

In Fig. |2] we show the chemical potential fi (in units of 
U) versus x = 2dt/U phase diagram for (a) two-dimensional 
and (b) three-dimensional hypercubic lattices, where in these 
figures we choose the interspecies interaction to be attractive 

V = -0.85U. Comparing Eqs. (jSj and Q with Eqs. O 
and ( fT2] i. we expect that the excited state of the system to 
be a paired superfluid for all 1/ < when t 0. This is 
clearly seen in the figure where the dotted lines correspond to 
phase boundary for the Mott insulator to superfluid state as de- 
termined from the third-order strong-coupling expansion, the 
hollow pink-squares correspond to the extrapolation fits for 
the single-particle and single-hole excitations (shown only for 
illustration purposes), and the solid black-circles correspond 
to the extrapolation fits for the two-particle and two-hole ex- 
citations (this is the expected transition) discussed in the text. 

At t = 0, the chemical potential width of all Mott lobes 
are V + U = 0.1 5J7, which is in contrast with the single- 
species BH model. As t increases from zero, the range of /i 
about which the ground state is a Mott insulator decreases here 
as well, and the Mott insulator phase disappears at a critical 
value of t, beyond which the system becomes a paired super- 
fluid. The strong-coupling expansion again overestimates the 
phase boundaries, and it again leads to unphysical pointed tips 
for all Mott lobes. In addition, a short list of V/U versus the 
critical points Xc = 2dtc/U are presented for the first two 
Mott insulator lobes in Table|T] Our results are consistent with 
the expectation that, for small V, the locations of the tips in- 
crease as a function of V, because the presence of a nonzero 

V is what allowed these states to form in the first place. How- 
ever, when V is larger than some critical value (~ 0.6U), the 
locations of the tips decrease, and they eventually vanish when 

V — —U. This may indicate an instability towards a collapse 
since at this point Uf-^Ui^ is exactly equal to U^^. 



Compared to the V > case shown in Fig. [T] note that 
shape of the Mott insulator to paired superfluid phase bound- 
ary is very different, showing a re-entrant behavior in all di- 
mensions from paired superfluid to Mott insulator and again 
to a paired superfluid phase, as a function of t. Our results 
are consistent with an early numerical time-evolving block 
decimation (TEBD) calculation yj], where such a re-entrant 
quantum phase transition in one dimension was predicted. 

The re-entrant quantum phase transition occurs when co- 
efficient of the hopping term in Eq. (fT2] l is negative [so 
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(a) Two dimensions (V=-0.85U) 
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(b) Three dimensions (V=-0.85U) 
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FIG. 2. (Color online) Chemical potential fi (in units of U) versus 
X = 2dt/U phase diagram for (a) two- and (b) three-dimensional 
hypercubic lattices with tf — — t, Uff — Un = U, Ufi = 
V = —0.8517, Tif = ni — n, and fj.^ — = p. The dotted lines 
correspond to phase boundary for the Mott insulator to superfluid 
state determined from the third-order strong-coupling expansion, the 
hollow pink-squares to the extrapolation fit for the single-particle or 
single-hole excitations (shown only for illustration purposes), and 
the solid black-circles to the extrapolation fit for the two-particle or 
two-hole excitations (the expected transition) discussed in the text. 



that the two-hole excitation branch has a negative slope in 
+ versus to- phase diagram when 0], i.e. 

C/fj.) + 2na{n„ + 1) /Uaa\zt^ term, which occurs for the first 
few Mott lobes beyond a critical U^^. When this coefficient 
is negative, its value is most negative for the first Mott lobe. 



TABLE II. List of the critical points (location of the tips) Xc = 
2dtc/U that are found from the chemical potential extrapolation 
technique described in the text. Here, =1^ = t, Uff = Un = U, 
Ufi = V, = ni = n, and /i-j- = /ij, = /i. These critical 
points for the two-particle or two-hole excitations are determined 
from Eqs. Jl lb and il2\ when < 0. Note that, for small V, Xc's 
tend to increase as a function of V, since the presence of a nonzero 

V is what allowed these states to form in the first place. However, 
Xc's decrease beyond a critical V, and they eventually vanish when 

V — —U, which may indicate an instability towards a collapse. 



V fU 


d = 2 


d = 3 


n = 1 


n = 2 


n = 1 


n = 2 


-0.01 


0.0543 


0.0337 


0.0611 


0.0379 


-0.03 


0.0937 


0.0582 


0.105 


0.0655 


-0.05 


0.121 


0.0749 


0.136 


0.0843 


-0.07 


0.142 


0.0883 


0.160 


0.0994 


-0.1 


0.169 


0.105 


0.190 


0.118 


-0.2 


0.233 


0.145 


0.262 


0.164 


-0.3 


0.277 


0.173 


0.311 


0.195 


-0.4 


0.307 


0.193 


0.345 


0.217 


-0.5 


0.325 


0.205 


0.366 


0.230 


-0.6 


0.331 


0.209 


0.372 


0.235 


-0.7 


0.321 


0.203 


0.362 


0.228 


-0.8 


0.291 


0.183 


0.327 


0.206 


-0.9 


0.225 


0.141 


0.253 


0.159 


-0.93 


0.193 


0.121 


0.217 


0.136 


-0.95 


0.166 


0.103 


0.187 


0.116 


-0.97 


0.1304 0.0812 


0.147 


0.0913 


-0.99 


0.0764 0.0474 


0.0860 0.0534 



and therefore the effect is strongest there. However, the coef- 
ficient increases and eventually becomes positive as a function 
of filling, and thus the re-entrant behavior becomes weaker as 
filling increases, and it eventually disappears beyond a critical 
filling. For the parameters used in Fig.|2] this occurs only for 
the first lobe, as can be seen in the figures. We also note that 
the sign of this coefficient is independent of the dimensional- 
ity of the lattice, since z — 2d enters into the coefficient only 
as an overall factor. 

What happens when t-^ ^ ti and/or U^^ ^ We donot 
expect any qualitative change for attractive interspecies inter- 
actions. However, for repulsive interspecies interactions, this 
lifts the degeneracy of the single-particle or single-hole exci- 
tation energies. While the transition is from a double Mott 
insulator to a double superfluid of both species in the degen- 
erate case, it is from a double-Mott insulator of both species 
to a Mott insulator of one species and a superfluid of the other 
in the nondegenerate case. 

V. CONCLUSIONS 

We analyzed the zero temperature phase diagram of the 
two-species Bose-Hubbard (BH) model with on-site boson- 
boson interactions in d-dimensional hypercubic lattices, in- 
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eluding both the repulsive and attractive interspecies in- 
teraction. We used the many-body version of Rayleigh- 
Schrodinger perturbation theory in the kinetic energy term 
with respect to the ground state of the system when the ki- 
netic energy term is absent, and calculate ground state ener- 
gies needed to carry out our analysis. This technique was 
previously used to discuss the phase diagram of the single- 
species BH model IH^IH, extended BH model iQ, and 
of the hardcore BH model with a superlattice [ZS*], and its 
results showed an excellent agreement with Monte Carlo sim- 
ulations ll23i I25I1 . Motivated by the success of this technique 
with these models, here we generalized it to the two-species 
BH model, hoping to develop an analytical approach which 
could be as accurate as the numerical ones. 

We derived analytical expressions for the phase boundary 
between the incompressible Mott insulator and the compress- 
ible superfluid phase up to third order in the hoppings. We also 
proposed a chemical potential extrapolation technique based 
on the scaling theory to extrapolate our third-order power se- 
ries expansion into a functional form that is appropriate for the 
Mott lobes. In particular, when the interspecies interaction is 
sufficiently large and attractive, we found a re-entrant quan- 
tum phase transition from paired superfluid (superfluidity of 
composite bosons, i.e. Bose-Bose pairs) to Mott insulator and 
again to a paired superfluid in all one, two and three dimen- 



sions. Since the available Monte Carlo calculations |'9','l^ do 
not provide the Mott insulator to superfluid transition phase 
boundary in the experimentally more relevant chemical po- 
tential versus hopping plane, we could not compare our results 
with them. This comparison is highly desirable to judge the 
accuracy of our strong-coupling expansion results. 

A possible direction to extend this work is to consider the 
limit where hopping of one-species is much larger than the 
other In this limit, the two-species BH model reduces to 
the Bose-Bose version of the Falicov-Kimball model [28], the 
Fermi-Fermi version of which has been widely discussed in 
the condensed-matter literature and the Fermi-Bose version 
has just been studied |29]. It is known for such models that 
there is a tendency towards both phase separation and density 
wave order |30], which requires a new calculation partially 
similar to that of Ref. ll24ll . One can also examine how the 
momentum distribution changes with the hopping in the insu- 
lating phases ll23l[3ll] . which has direct relevance to ultracold 
atomic experiments. 
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